function [delta_ig, a_star_ig, z_star_ig, a_ig, aM_g, AMs, ADs, betaT, beta_s, alphaL_s, alphaI_s, alpha_ks, barL_rs, Z_rs, D, E_s, F] = ... 
    invert_param_full(x_ig, m_ig, tau_star_ig, tau_ig, trade_con_share, tot_labor_comp, tot_intermediate, sales_s, IO_sales, Lsr, ...
    Dsg, DX, DM, omega_star, sigma, eta, kappa, sigma_star, gammaM, gammaX, tol)

%This script implements the inversion of model parameters using observed variables
%Corresponding author: Rodrigo Adao
%Date: 09/11/2024

[N, G] = size(x_ig);
[R, S] = size( Lsr );

%Variety-level export parameters
delta_ig = ones(size(x_ig));
a_star_ig = ( (1 + tau_star_ig).^(sigma_star) ).*x_ig;

%Variety-level import parameters
z_star_ig = ( ( (1 + tau_ig).*m_ig ).^(-omega_star) )./(1 + tau_ig);

a_ig =   (1+ tau_ig).*m_ig ;
a_ig_aux = sum( a_ig , 1) ; 
a_ig(:,a_ig_aux<tol) = 1;
a_ig_aux(a_ig_aux<tol) = N;
a_ig = a_ig*spdiags(1./a_ig_aux', 0, G, G);
if max( abs(sum(a_ig,1) - 1) ) > tol 
    error('sum of a_ig across i is not one for some g')
end

varphi_ig = ( a_ig.^(1+gammaM*omega_star) ).*( z_star_ig.^(gammaM*(1-sigma)) ).*( (1 + tau_ig).^(gammaM*(1-sigma)) );
varphi_ig = varphi_ig.^( 1/(1+omega_star*sigma*gammaM) ) ;

%Product-level import parameters
m_g = sum( (1+tau_ig).*m_ig , 1)';
varphi_g = sum( varphi_ig , 1)';

aM_g = ( m_g.^( (1+omega_star*eta*gammaM)/(1+omega_star*gammaM) ) ).*( varphi_g.^( - ( (1-eta)*(1+omega_star*sigma*gammaM) )/( (1-sigma)*(1+omega_star*gammaM) ) ) );
aM_g(m_g<tol) = 0;   %assign parameter = 0 if import of product g = 0
aM_gs_aux = ( ( Dsg*aM_g )'*Dsg )'; 
aM_g(aM_gs_aux<tol) = 1; %wlog, assign equal shifters to all products if there are no imports in a sector
aM_gs = ( ( Dsg*aM_g )'*Dsg )';

%adjust size 
aM_g = aM_g./aM_gs;
if max( abs(Dsg*aM_g - 1) ) > tol 
    error('sum of aM_g across g is not one for some s')
end

%sector-level production parameters
alphaL_s = tot_labor_comp./sales_s;
alphaI_s = tot_intermediate./sales_s;
alphaK_s = 1 - alphaL_s - alphaI_s; 
alpha_ks = IO_sales*diag(1./tot_intermediate);
if min( alphaK_s ) < 0
    error('alphaKs negative for some sector')
end
if max( abs(sum(alpha_ks,1) - 1) ) > tol 
    error('sum of alpha_ks across k is not one for some sector')
end

%Sector-level final demand parameters
EXs = Dsg*sum(x_ig,1)';
Ms = Dsg*m_g;
F_s = Ms - EXs + ( eye(S) - alpha_ks*diag(alphaI_s) )*sales_s ;
%F_s = ( - EXs + (eye(S) - diag(ADs)*alpha_ks*diag(alphaI_s) )*sales_s )./ADs;
if min(F_s) < 0 
    error('final spending is negative for some sector')
end

betaT = sum(trade_con_share);
betaNT = 1 - betaT;
beta_s = F_s/sum(F_s);
if abs(sum(beta_s) - 1) > tol 
    error('sum of beta_s across s is not one')
end
%corr([beta_s trade_con_share/betaT])

%Compute spending by sector
E_s = F_s + ( alpha_ks*diag(alphaI_s) )*sales_s;

%Sector-level import parameters
AMs = (Dsg*m_g)./E_s ;
ADs = 1 - AMs;
if max( abs(AMs + ADs - 1)) > tol 
    error('AM_s + AD_s is not one for some s')
end

%Check with sales and spending are consistent
Y_s = ( inv(eye(S) - diag(ADs)*alpha_ks*diag(alphaI_s) ) )*( ADs.*F_s + EXs );
if max(abs(Y_s./sales_s - 1)) > tol
    error('gross sales do not match demand')
end
W_s = Y_s.*alphaL_s;

%sector-region parameters
barL_rs = Lsr*spdiags( W_s./sum(Lsr,1)' , 0, S, S );
if max(abs(sum(barL_rs,1)./W_s' -1)) > tol 
    error('sum of barL_sr across r is not equal to aggregate labor payments for some sector')
end

Z_rs = barL_rs*spdiags( 1./alphaL_s, 0, S, S ) ;
Z_rs = Z_rs.^repmat(alphaK_s', R, 1);

%International Transfer
RevenueM = DM*sum( tau_ig.*m_ig, 'all' ); 
IncomeT = (1-alphaI_s)'*Y_s;
RevenueX = DX*sum( tau_star_ig.*x_ig, 'all' );
F = sum(F_s)/betaT;
D = (1-betaNT)*F - IncomeT - RevenueM - RevenueX; 

%Final checks
income = (D + RevenueM + RevenueX + IncomeT)/(1-betaNT);

%sector spending
Es_check = beta_s*betaT*income + alpha_ks*diag(alphaI_s)*Y_s;
if max( abs(Es_check./E_s - 1)) > tol 
    error('sector spending do not match demand')
end

%sector output
Ys_check = ADs.*E_s + Dsg*sum(x_ig,1)';
if max( abs(Ys_check./Y_s - 1)) > tol 
    error('sector output do not match demand')
end

%sector imports
 Imp_check = Dsg*m_g;
 if max( abs( (AMs.*E_s)./Imp_check -1 ) ) > tol 
     error('sector imports do not match demand')
 end

%sector market clearing
supply_s = (Y_s - Dsg*sum(x_ig,1)' + Dsg*sum( (1+tau_ig).*m_ig,1)' );
if max( abs( supply_s./E_s -1)) > tol 
    error('sector market clearing does not hold')
end


end

